home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / BSSTEP.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  60 lines
  1. PROCEDURE bsstep(VAR y: glyarray; dydx: glyarray; nv: integer; VAR x: real;
  2.       htry,eps: real; yscal: glyarray; VAR hdid,hnext: real);
  3. (* Programs using routine BSSTEP must define the type
  4. TYPE
  5.    glyarray = ARRAY [1..nv] OF real;
  6. in the main routine. *)
  7. LABEL 99;
  8. CONST
  9.    imax=11;
  10.    nuse=7;
  11.    one=1.0e0;
  12.    shrink=0.95e0;
  13.    grow=1.2e0;
  14. VAR
  15.    j,i: integer;
  16.    xsav,xest,h,errmax: real;
  17.    ysav,dysav,yseq,yerr: glyarray;
  18.    nseq: ARRAY [1..imax] OF integer;
  19. BEGIN
  20.    nseq[1] := 2; nseq[2] := 4; nseq[3] := 6;
  21.    nseq[4] := 8; nseq[5] := 12; nseq[6] := 16;
  22.    nseq[7] := 24; nseq[8] := 32; nseq[9] := 48;
  23.    nseq[10] := 64; nseq[11] := 96;
  24.    h := htry;
  25.    xsav := x;
  26.    FOR i := 1 TO nv DO BEGIN
  27.       ysav[i] := y[i];
  28.       dysav[i] := dydx[i]
  29.    END;
  30.    WHILE true DO BEGIN
  31.       FOR i := 1 TO imax DO BEGIN
  32.          mmid(ysav,dysav,nv,xsav,h,nseq[i],yseq);
  33.          xest := sqr(h/nseq[i]);
  34.          rzextr(i,xest,yseq,y,yerr,nv,nuse);
  35.          errmax := 0.0;
  36.          FOR j := 1 TO nv DO BEGIN
  37.             IF (errmax < abs(yerr[j]/yscal[j])) THEN
  38.                errmax := abs(yerr[j]/yscal[j])
  39.          END;
  40.          errmax := errmax/eps;
  41.          IF (errmax < one) THEN BEGIN
  42.             x := x+h;
  43.             hdid := h;
  44.             IF (i = nuse) THEN hnext := h*shrink
  45.             ELSE IF (i = nuse-1) THEN hnext := h*grow
  46.             ELSE hnext := (h*nseq[nuse-1])/nseq[i];
  47.             GOTO 99
  48.          END
  49.       END;
  50.       h := 0.25*h;
  51.       IF (((imax-nuse) DIV 2) > 0) THEN BEGIN
  52.          FOR i := 1 TO ((imax-nuse) DIV 2) DO h := h/2
  53.       END;
  54.       IF ((x+h) = x) THEN BEGIN
  55.          writeln('pause in routine BSSTEP');
  56.          writeln('step size underflow'); readln
  57.       END
  58.    END;
  59. 99:   END;
  60.